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ABSTRACT 

Numerical simulations of Milky- Way size Cold Dark Matter (CDM) halos predict a steeply rising 
mass function of small dark matter subhalos and a substructure count that greatly outnumbers the 
observed satellites of the Milky Way. Several proposed explanations exist, but detailed comparison 
between theory and observation in terms of the maximum circular velocity (Vmax) of the subhalos is 
hampered by the fact that Vmax for satellite halos is poorly constrained. We present comprehensive 
mass models for the well-known Milky Way dwarf satellites, and derive likelihood functions to show 
that their masses within 0.6 kpc (Mo. 6) are strongly constrained by the present data. We show that 
the Mo. 6 mass function of luminous satellite halos is flat between ~ 10 7 and 10 8 Af Q . We use the 
"Via Lactea" N-body simulation to show that the Ma. 6 mass function of CDM subhalos is steeply 
rising over this range. We rule out the hypothesis that the 11 well-known satellites of the Milky Way 
are hosted by the 11 most massive subhalos. We show that models where the brightest satellites 
correspond to the earliest forming subhalos or the most massive accreted objects both reproduce the 
observed mass function. A similar analysis with the newly-discovered dwarf satellites will further test 
these scenarios and provide powerful constraints on the CDM small-scale power spectrum and warm 
dark matter models. 

Subject headings: Cosmology: dark matter, theory-galaxies 



1. INTRODUCTION 

It is now well-established that numerical simulations 
of Cold Dark Matter (CDM) halos predict many orders 
of magnitude more dark matter subhalos around Milky 
Wa y-sized galaxies than observed dwarf satellite galax- 
ies (IKlvpin et alJl999HMoore et al.ll999trDiemand et all 
12007( 1 . Within the context of the CDM paradigm, there 
are well-motivated astrophysical solutions to this 'Miss- 
ing Satellites Problem' (MSP) that rely on the idea 
that galaxy fo rmation is inefficient in the smallest dark 
matter halos (IBullock et all 120001: iBenson et all | 200 " 
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Somer ville 2002; St oehr et al.ll2002l:lKravtsov et al.l l200l 
Moore et al.ll2006t iGnedin fc Kravtsovl [20061 ). However, 
from an observational perspective, it has not been possi- 
ble to distinguish between these solutions. 

A detailed understanding of the MSP is limited by our 
lack of a precise ability to characterize the dark matter 
halos of satellite galaxies. From an observational per- 
spective, the primary constraints come from the veloc- 
ity dispersion of ~ 200 stars in the population of dark 
matter dominated dwarf spheroidal gal a xies (dSphs) 
( Wilk inson et~ail l2004t lLokas et all 120051: iMunoz et al l 
2005112006: Westfall et al. 2006; Walker et al.ll2005ll2006f 



Sohn et al.ir2006t iGilmore et al.ll2007h . However, the ob- 



served extent of the stellar populations in dSphs are ~ 
kpc, so these velocity dispersion measurements are only 
able to probe properties of the halos in this limited radial 
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regime. 

From the theoretical perspective, dissipationless nu- 
merical simulations typically characterize subhalo counts 
as a function of the total bound mass or maximum cir- 
cular velocity, Vmax- While robustly determined in sim- 
ulations, global quantities like Vmax are difficult to con- 
strain observationally because dark halos can extend well 
beyond the stellar radius of a satellite. Indeed stellar 
kinematics alone provide only a lower limit on the halo 
Vnax value (see below). This is a fundamental limita- 
tion of stellar kinematics that cannot be remedied by 
increasing the n umber of stars used in the velocity dis- 
persion analysis ([Strigari et al.| [2007). Thus determining 
V m ax values for satellite halos requires a theoretical ex- 
trapolation. Any extrapolation of this kind is sensitive 
to the predicted density structure of subhalos, which de- 
pends on cosmology, power -spectrum normalization, and 
the nature of dark matter ([Zentner &; Bullockll2003l ). 

Our inability to determinate Vnax is the primary lim- 
itation to test solutions to the MSP. One particular so- 
lution states that the masses of the dSphs have been 
systematically underestimated, so that the ~ 10 bright- 
est satellites r eside systema ticall y in the ~ 10 m ost mas- 
sive subhalos ([Stoehr et al.ll2002tlHavashi et alj|2003) . A 
byproduct of this solution is that there must be a sharp 
mass cutoff at some current subhalo mass, below which 
galaxy formation is suppressed. Other solutions, based 
on reionization suppression, or a characteristic halo mass 
scale prior to subhalo accretion predict that the sup- 
pression comes in gradually with cu rrent subhalo mass 
(IBullock et all 120001 : iKravtsov etai] 12004 iMoore et all 
I2006D . 

In this paper, we provide a systematic investigation of 
the masses of the Milky Way satellites. We highlight 
that in all cases the total halo masses and maximum 
circular velocities are not well-determined by the data. 
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We instead use the fact that the total cumulative mass 
within a characteristic radius ~ .6 kpc is much bet - 
ter determined by the present data (|Strigari et alJ feOO?) . 
We propose using this mass, which we define as M .6, 
as the favored means to compare to the dark halo pop- 
ulation in numerical simulations. Unlike U max , Mo. 6 is 
measured directly and requires no cosmology-dependent 
or dark-matter- model-dependent theoretical prior. 

In the following sections, we determine the M 6 mass 
function for the Milky Way satellites, and compare it 
to the corresponding mass function measured directly in 
th e high-resolution "Via L actea" substructure simulation 
of iDiemand et al.1 (|2007f L We rule out the possibility 
that there is a one-to- one correspondence between the 11 
most luminous satellites and the most massive subhalos 
in Via Lactea. We find that MSP solutions based on 
rcionization and/or characteristic halo mass scales prior 
to accretion are still viable. 

2. MILKY WAY SATELLITES 

Approximately 20 satellite galaxies can be classified 
as residing in MW subhalos. Of these, ~ 9 were 
discovered within the last two years and have very 
low luminosities and surface brightnesses (IWillman et al.l 
20051: IWillman et all 120051: iBelokurov et al.ll2006l l2007t 
Zucker et al. 2006allj _ The lack of precision in these 
numbers reflects the ambiguity in the classification of the 
newly-discovered objects. The nine 'older' galaxies clas- 
sified as dSphs are supported by thei r velocity disp ersion, 
and exhibit no rotational motion (iMateol Il998f ). Two 
satellite galaxies, the Small Megallanic Cloud (SMC) 
and Large Megallanic Cloud (LMC) , are most likely sup- 
ported by some combination of rotation and dispersion. 
Stellar kinematics suggest that the LMC and SMC are 
likely the most massive satellite systems of the Milky 
Way. 

We focus on determining the masses of the nine most 
well-studied dSphs. The dark matter masses of the dSphs 
are determined from the line-of-sight velocities of the 
stars, which trace the total gravitational potential. We 
assume a negligible contribution of the stars to the grav- 
itational potential, which we find to be an excellent ap- 
proximation. The dSph with the smallest mass-to-light 
ratio is Fornax, though even for this case we find that the 
star s generally do not affe ct the dynamics of the system 
see lLokail2002l: IWul[2007l and below). 

Under the assumptions of equilibrium and spherical 
symmetry, the radial component of the stellar velocity 
dispersion, ay, is linked to the total gravitational poten- 
tial of the system via the Jeans equation, 



spherically-symmetric King profile (|Kindll962l ). 



dr 



(1) 



Here, v± is the stellar density profile, V 2 = GM{r)/r 
includes the total gravitating mass of the system, and 
the parameter (3{r) = \ — a 2 /a 2 characterizes the dif- 
ference between the radial (oy) and tangential (ag) ve- 
locity dispersions. The observable velocity dispersion is 
constructed by integrating the three- dimensional stellar 
radial velocity dispersion profile along the line-of-sight, 
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where R is the projected radius on the sky. The surface 
density of stars in all dSphs are reasonably well-fit by a 
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where r t and Thing are fitting parameters denoted as the 
tidal and core radii. 6 The spherically symmetric stellar 
density can be obtained with an integral transformation 
of the surface density. From the form of Eq. ([2]), the 
normalization in Eq. ([3]) is irrelevant. 

Some dSphs show evidence for multiple, dynamically 
distinct stellar populations, with each population de- 
scribed by its own surface density and velocity dispersion 
(iHarbeck et al.ll2001fc iTolstov et all 12004 IWestfajl et all 
120061 : iMcConnachie et al.l 120061 : llbata et al.ll2006l ). In a 
dSph with i = l...N p populations of stars, standard ob- 
servational methods will sample a density-weighted av- 
erage of the populations: 
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where vi and Oi are the density profile and radial stellar 
velocity dispersion of stellar component i. In principle, 
each component has its own stellar velocity anisotropy 
profile, fli{r). In this case, Equation ([1} is valid for v± and 
ay as long as an effective velocity anisotropy is defined 
as 



From these definitions, we also have ^ z /* t a 
ho 2 
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The conclusion we draw from this argument is that 
the presence of multiple populations will not affect the 
inferred mass structure of the system, provided that the 
velocity anisotropy is modeled as a free function of ra- 
dius. Since the functional form of the stellar velocity 
anisotropy is unknown, we allow for a general, three pa- 
rameter model of the velocity anisotropy profile, 



(7) 



A profile of this form allows for the possibility for (3{r) to 
change from radial to tangential orbits within the halo, 
and a constant velocity anisotropy is recovered in the 
limit Poo — > 0, and (3q — > const. 

In Equations |T]) and ([2]), the radial stellar velocity dis- 
persion, oy, depends on the total mass distribution, and 
thus the parameters describing the dark matter density 
profile. Dissipation-less N-body simulations show that 
the density profiles of CDM halos can be characterized 
as 

7-T(l + r)° 

where r s and p s set a radial scale and density normaliza- 
tion and 7 and 8 parameterize the inner and outer slopes 

6 Our results are insensitive to this particular parameterization 
of the light profile. 
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of the distribution. For dark matter halos unaffected 
by tidal interactions, the most recent high-resolution 
simulations find 5 + 7 « 3 works well for the outer 
slope, while 0.7 < 7 < 1.2 works well down to ~ 0.1% 
of ha lo virial radii () Navarro et al.l 120041 : iDiemand et al.l 
2005). This interior slope is not altered as a subhalo 
loses mass from tidal inte ractions with the MW potential 
(|Kazantzidis et all I2004T ). The outer slope, 5, depends 
more sensitively on the tidal interactions in the halo. The 
majority of the stripped material will be from the outer 
parts of halos, and thus S of subhalo density profiles will 
become steeper than those of field halos. Subhalos are 
characterized by outer slopes in the range 2 < 6 < 3. 

Given the uncertainty in the j3(r) and p(r) profiles, 
we are left with nine shape parameters that must be 
constrained via line-of-sight velocity dispersion measure- 
ments: p s , r s , 0o, Poo, rp, 7, S, rking, and r t . While the 
problem as posed may seem impossible, there are a num- 
ber of physical parameters, which are degenerate between 
different profile shapes, that are well constrained. Specif- 
ically, the stellar populations constrain V c (r) within a 
radius comparable to the stellar radius r t ~ kpc. As 
a result, quantities like the local density and integrated 
mass with in the stellar radius are determined with high 
precision (jStrigari et al.|[2006h . while quantities that de- 
pend on the mass distribution at larger radii are virtually 
unconstrained by the data. 

It is useful to determine the value of the radius where 
the constraints are maximized. The location of this char- 
acteristic radius is determined by the form of the integral 
in Eq. (J5J). We can gain some insight using the exam- 
ple of a power-law stellar distribution f*(r), power-law 
dark matter density profile p oc r -7 * , and constant ve- 
locity anisotropy. The line-of-sight velocity dispersion 
depends on the three-dimensional stellar velocity disper- 
sion, which can be written as 



Gv*{r)M{r)r 2f} - 2 dr cx r 2 ^* (9) 



From the shape of the King profile, the majority of 
stars reside at projected radii r^ng < R < r t , where 
the stellar distribution is falling rapidly v+ ~ r -3.5_ 
In this case, for (3 = 0, the LOS component scales as 
a L( R ) <* Ir r~ - 5 -^ (r 2 - i? 2 )" 1 / 2 dr and is dominated 
by the mass and density profile at the smallest relevant 
radii, r ~ r^ng- For R < rki ng , v* oc r^ 1 and af os is simi- 
larly dominated by r ~ rki ng contributions. We note that 
although the scaling arguments above hold for constant 
velocity anisotropies, they can be extended by consider- 
ing anisotropies that vary significantly in radius. They 
are also independent of the precise form of I*, provided 
there is a scale sim i lar to rking- 

In IStrigari et alj (|2007f) , it was shown that typical ve- 
locity dispersion profiles best constrain the mass (and 
density) within a characteristic radius ~ 2 rking ■ For 
example, the total mass within 2rking is constrained to 
within ~ 20% for dSphs with ~ 200 line-of-sight veloc- 
ities. Note that when deriving constraints using only 
the innermost stellar velocity dispersion and fixing the 
anisotropy to (3 — 0, the characterist ic radius for best 
const raints decreases to ~ 0-5rki ng (e.g. lPenarrubia et al.1 
12001 . 

As listed in Table [TJ the Milky Way dSphs are ob- 
served to have variety of rking values, but rking ~ 0.3 



kpc is t ypical . The values of rking and r t are taken from 
iMateol (|1998f ). In order to facilitate comparison with 
simulated subhalos, we chose a single characteristic ra- 
dius of 0.6 kpc for all the dwarfs, and we represent the 
mass within this radius as Mq.q = M(< 0.6 kpc). The 
relative errors on the derived masses are unaffected for 
small variations in the characteristic radius in the range 
~ 1.5 — 2.5 rking- Deviations from a true King profile at 
large radius (near rt) do not affect these arguments, as 
long as there is a characteristic scale similar to rking in 
the surface density. The only dSph significantly affected 
by the choice of 0.6 kpc as the characteristic radius is 
Leo II, which has r t — 0.5 kpc. Since the characteristic 
radius is greater than twice rking-, the constraints on its 
mass will be weakest of all galaxies (with the exception 
of Sagittarius, as discussed below). 

3. DARK MATTER HALO MASSES AT THE 
CHARACTERISTIC RADIUS 



We use the following data s ets: [W ilkinson 'et al.1 
d200l:lMunoz et al.l (127)051 [20061) : IWestfall et all rt2006h: 
iWalker et all (120051 . 120061) : iSohn et all ( |2006h : Sieeal et 
al. in preparation. These velocity dispersions are deter- 
mined from the line-of-sight velocities of ~ 200 stars in 
each galaxy, although observations in the coming years 
will likely increase this number by a factor ~ 5—10. 
From the data, we calculate the x 2 , defined in our case 
as 



X 



1 

\&obs,i 

t=l 



(10) 



Here o 2 ohs is the observed velocity dispersion in each bin, 
<j 2 h is the theoretical value, obtained from Eq. and 
e 2 are errors as determined from the observations. 

It is easy to see that, when fitting to a single data set 
of ~ 200 stars, parameter degeneracies will be signifi- 
cant. However, from the discussion in section [21 Mo. 6 
is well-determined by the LOS data. To determine how 
well M .6 is constrained, we construct likelihood func- 
tions for each galaxy. When thought of as a function 
of the theoretical parameters, the likelihood function, C, 
is defined as the probability that a data set is acquired 
given a set of theoretical parameters. In our case £ is a 
function of the parameters 7, 6, r s , p s , and /3q, /3oo, rp, 

and is defined as L — e~ x I 2 . In writing this likelihood 
function, we assume that the errors on the measured ve- 
locity dispersions are Gaussian, which we find to be an 
excellent approxima tion to the errors for a dSph system 
(|Strigari et al.ll200l . We marginalize over the parame- 
ters 7, 5, r s , p s , and 0q, fiao, rp at fixed Afo.67 and the 
optimal values for Mq q are determined by the maximum 
of C . 

We determine C for all nine dSphs with velocity dis- 
persion measurements. For all galaxies we use the full 
published velocity dispersion profiles. The only galaxy 
that does not have a published velocity dispersion profile 
is Sagittarius, and for this galaxy w e use the central ve- 
locity dispersion from lMateol (| 19981 ) . The mass modeling 
of Sagittarius is further complicated by the f act that it is 
experiencing tidal interactio ns with the MW (|Ibata et al.l 
1997; Majewski ct al. 2003), so a mass estimate from the 
Jeans equation is not necessarily reliable. We caution 
that in this case the mass we determine is likely only an 
approximation to the total mass of the system. 
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Fig. 1. — The likelihood functions for the mass within 0.6 kpc 
for the nine dSphs, normalized to unity at the peak. 
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Fig. 2. — The velocity dispersion for Ursa Minor as a function of 
radial distance, along with the model that maximizes the likelihood 
function. 



We determine the likelihoods by marginalizing over 
the following ranges of the velocity anisotropy, inner 
and outer slopes: —10 < (3q < 1, —10 < (3^ < 1, 
0.1 < r < 10 kpc, 0.7 < 7 < 1.2, and 2 < S < 3. 
As discussed above, these ranges for the asymptotic in- 
ner and outer slopes are appropriate because we are con- 
sidering CDM halos. It is important to emphasize that 
these ranges are theoretically motivated and that obser- 
vations alone do not demand such restrictive choices. 
It is possible to fit all of the dSphs at present with 
a constant density cores with scale- leng ths ~ 100 pc 
(jStrigari et al.ll2006t iGilmore et all 120071 ). although the 
data by no means demand such a situation. Though 
we consider inner and outer slopes in the ranges quoted 
above, our results are not strongly affected if we widen 
these intervals. For example, we find that if we allow 
the inner slope to vary down to 7 = 0, the widths of the 
likelihoods are only changed by ~ 10%. This reflects the 
fact that there is a negligible degeneracy between Mo. 6 
and the inner and outer slopes. 

We are left to determine the regions of p s and r s pa- 
rameter space to marginalize over. In all dSphs, there 
is a degeneracy in this parameter space, telling us that 
it is not possible to derive an upp er limit on this pai r 
of parameters from the data alone (|Strigari et al.l f2006). 
While this degeneracy is not important when determin- 
ing constraints on Mo. 6, it is the primary obstacle in 
determining V mayi . From the fits we present below, we 
find that the lowest r s value that provides an acceptable 
fit is ~ 0.1 kpc, and we use this as the lower limit in all 
cases. In our fiducial mass models, we conservatively re- 
strict the maximum value of r s using the known distance 
to each dSph. In this case, we use 0.1 kpc < r s < D/2, 
where D is the distance to the dSph. 

In Figure Q] we show the Mq.q likelihood functions for 
all of the dSphs. As is shown, we obtain strong con- 
straints on Mo. 6 in all cases except Sagittarius, for which 
we use only a central velocity dispersion. Table [1] sum- 
marizes the best fitting Mo. 6 values for each dwarf. The 
quoted errors correspond to the points where the likeli- 
hood falls to 10% of its peak value. The upper panel of 



Figure [3] shows Mo. 6 values for each dwarf as a function 
of luminosity. In Figure O we show an example of the 
velocity dispersion data as a function of radial distance 
for Ursa Minor, along with the model that maximizes 
the likelihood function. For all galaxies, we find \ 2 P er 
degree of freedom values < 1 . 

The maximum likelihood method also allows us to con- 
strain the mass at other radii spanned by the stellar dis- 
tribution. The sixth column of Table 1 provides the inte- 
grated mass within each dwarf's King tidal radius. This 
radius roughly corresponds to the largest radius where 
a reasonable mass constraint is possible. As expected, 
the mass within r t is not as well determined as the mass 
within 2 r^i ng . From these masses we are able to deter- 
mine the mass-to- light ratios within r t , which we present 
in the seventh column of Table [TJ In the bottom panel 
of Figure [3l we show mass-to-light ratios within r t as a 
function of dwarf luminosity. We see the standard result 
that the observable mass-to-light ratio increases with de- 
creasing luminosity (|Mateol[l998h . Note, however, that 
our results are inconsistent with the idea that all of the 
dwarfs have the same integrated mass within their stellar 
extent. We note that for Sagittarius, we can only obtain 
a lower limit on the total mass-to-light ratio. 

The last two columns in Table 1 list constraints on 
Knax for the dSphs. Column 8 shows results for an anal- 
ysis with limits on r s as described above. In this case, the 
integrated mass within the stellar radius is constrained 
by the velocity dispersion data, but the halo rotation ve- 
locity curve, V c (r), can continue to rise as r increases 
beyond the stellar radius in an unconstrained manner. 
The result is that the velocity dispersion data alone pro- 
vide only a lower limit on V^ lax . 

Stronger constraints on y ma x can be obtained if we 
limit the range of r s by imposing a cosmology-dependent 
prior on the dark matter mass profile. CDM simula- 
tions have shown that there is a correlation between 
Knax and r max for halos, where r max is the radius where 
the circular velocity peaks. Because subhalo densities 
will depend on the collapse time and orbital evolution 
of each system, the precise Vmax-rmax relation is sen- 
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TABLE 1 

Parameters Describing Milky Way Satellites. 
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0.20 


0.80 
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0.39 
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44+31 
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25+5 
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0.43 


3 4+0-7 
l.Q 


e 7+2.3 
D -'-2.5 
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68 
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14+2 


Sextans 
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4.01 
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Note. — Determination of the mass within 0.6 kpc and the maximum circular velocity for the dark matter halos of the dSphs. The errors are 
determined as the location where the likelihood function falls off by 90% from its peak value. For Sagittarius, no reliable estimate of V ma x with the 
CDM prior could be determined. The CDM prior is determined using the concordance cosmology with erg — 0.74, n — 0.95 (see text for details). 

values for the assumed prior. Note that in most cases, 
this prior degrades the quality of the fit, and the likeli- 
hood functions peak at a lower overall value. The mag- 
nitude of this effect is not large except for the cases of 
Leo II and Sagittarius. For Leo II, the peak likelihood 
with the prior occurs at a value that is below the 10% 
likelihood for the case without a prior on r s (i.e. the 
data seem to prefer a puffier subhalo than would be ex- 
pected in CDM). For Sagittarius, we are unable to obtain 
a reasonable fit within a subhalo that is typical of those 
expected. This is not too surprising. Sagittarius is be- 
ing tidally disrupted and its dark matter halo is likely 
atypical. 

We emphasize that the Vmax determinations listed in 
Column 9 are driven by theoretical assumptions, and can 
only be fairly compared to predictions for this specific 
cosmology (LCDM, cr 8 = 0.74). The M Q . 6 values in Col- 
umn 5 are applicable for any theoretical model, including 
non-CDM models, or CDM models with any normaliza- 
tion or power spectrum shape. 
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Fig. 3. — The mass within 0.6 kpc (upper) and the mass-to- 
light ratios within the King tidal radius (lower) for the Milky Way 
dSphs as a function of dwarf luminosity. The error-bars here are 
defined as the locations where the likelihoods fall to 40% of the 
peak values (corresponding to ~ la errors). The lines denote, 
from top to bottom, constant values of mass of 10 7 ,10 8 ,10 9 M Q . 

sitive to cosmology (e.g. cr s ) and the formation his- 
tory of the host halo itself (e.g. [Zentner fc Bullockll2003t 
Powerl2003l; iKa zantzidis et al. 2 0O4uBulk>ck fc Johnston] 
2005HBullock fc Johnstonll2006h . When converted to the 
relevant halo parameters, the imposed V max -r max rela- 
tion can be seen as a theoretical prior on CDM ha- 
los, restricting the parameter space we need to inte- 
grate over. In order to illustrate the technique, we adopt 
logio(^max) = 1.35(log 10 (Kiax/kms _1 ) - 1) - 0.196 kpc 
with a scatter of 0.2 in log 10 , as measured from simulated 
subh alos within the Via Lactea host halo (jDiemand et al.l 
2007). This simulation is for a LCDM cosmology with 
erg = 0.74 and n — 0.95. The scatter in the subhalo mass 
function increases at the very high mass end, which re- 
flects the fact that these most massive subhalos are those 
that are accreted most rece ntly (jZentner fc Bullockll2003t 
Ivan den Bosch et al.l l2005h . However, as we show below 
our results are not strongly dependent on the large scat- 
ter at the high mass end. 

Column 9 in Table 1 shows the allowed subhalo Vmax 



4. COMPARISON TO NUMERICAL SIMULATIONS 

The recently-completed Via Lactea run is the highest- 
resolution simulation of galactic substructure to date, 
containing an or der of magnitude mo re particles than 
its predecessors (|Diemand et al.l l2007| ) . As mentioned 
above, Via Lactea assumes a LCDM cosmology with 
o~ 3 = 0.74 and n = 0.95. For a d e tailed description of 
the simulation, see iDiemand et all (j2007). For our pur- 
poses, the most important aspect of Via Lactea is its 
ability to resolve the mass of subhalos on length scales 
of the characteristic radius 0.6 kpc. In Via Lactea, the 
force resolution is 90 pc and the smallest well-resolved 
length scale is 300 pc, so that the mass within 0.6 kpc 
is well-resolved in nearly all subhalos. Due to the choice 
of time steps we expect the simulation to underestimate 
local densities in the densest regions (by about 10% at 
densities of 9 x 10 7 M©/kpc 3 ). There is only one sub- 
halo with a higher local density than this at 0.6 kpc. For 
this subhalo, p(r = 0.6 kpc) = 1.4 x 10 8 M o /kpc 3 , so 
its local density might be underestimated by up to 10%, 
and the errors in the enclosed mass might be ~ 20% 
(jDiemand et al.ll2005| ). For all other subhalos the densi- 
ties at 0.6 kpc are well below the affected densities, and 
the enclosed mass should not be affected by more than 
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Fig. 4. — The mass within 0.6 kpc versus the maximum circular 
velocity for the mass ranges of Via Lactea subhalos corresponding 
to the population of satellites we study. 

10% by the finite numerical resolution. 

We define subhalos in Via Lactea to be the self-bound 
halos that lie within the radius R200 = 389 kpc, where 
R200 is defined to enclose an average density 200 times 
the mean matter density. We note that in comparing 
to the observed MW dwarf population, we could have 
conservatively chosen subhalos that are restricted to lie 
within the same radius as the most distant MW dSph 
(250 kpc). We find that this choice has a negligible effect 
on our conclusions - it reduces the count of small halos 
by ~ 10%. 

In Figure [3J we show how Ma. 6 relates to the more 
familiar quantity V max in Via Lactea subhalos. We note 
that the relationship between subhalo M 6 and V max will 
be sensitive to the power spectrum shape a nd normaliza- 
tion, as well as the nature of d ark matter (|Bullock et al.l 
120011: IZentner fc Bullockl I2003D . The relationship shown 
is only valid for the Via Lactea cosmology, but serves as 
a useful reference for this comparison. 

Given likelihood functions for the dSph M .6 values, 
we are now in position to determine the Mo. 6 mass func- 
tion for Milky Way (MW) satellites and compare this to 
the corresponding mass function in Via Lactea. For both 
the observations and the simulation, we count the num- 
ber of systems in four mass bins from 4 x 10 6 < Mo. 6 < 
4 x 10 8 M . This mass range is chosen to span the Mo. 6 
values allowed by the likelihood functions for the MW 
satellites. We assume that the two non-dSph satellites, 
the LMC and SMC, belong in t he highest mass bin, cor- 
responding to Mq.o > 10 s M (|Harris fc Zaritskvl [20061 : 
Ivan der Marel et alj|2003 ). 

In Figure [5] we show resulting mass functions for MW 
satellites (solid) and for Via Lactea subhalos (dashed, 
with Poisson error-bars). For the MW satellites, the cen- 
tral values correspond to the median number of galaxies 
per bin, which are obtained from the maximum values 
of the respective likelihood functions. The error-bars 
on the satellite points are set by the upper and lower 
configurations that occur with a probability of > 10~ 3 
after drawing 1000 realizations from the respective like- 
lihood functions. As seen in Figure[5l the predicted dark 
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Fig. 5. — The Mo. 6 mass function of Milky Way satellites and 
dark subhalos in the Via Lactea simulation. The red (short-dashed) 
curve is the total subhalo mass function from the simulation. The 
black (solid) curve is the median of the observed satellite mass 
function. The error-bars on the observed mass function represent 
the upper and lower limits on the number of configurations that 
occur with a probability of > 10 — 3 . 

subhalo mass function rises as ~ Mqq while the visi- 
ble MW satellite mass function is relatively flat. The 
lowest mass bin (Mo. 6 ~ 9 x 1O 6 M0) always contains 1 
visible galaxy (Sextans). The second-to- lowest mass bin 
(Ma. 6 ~ 2.5 x lO 7 Af ) contains between 2 and 4 satellites 
(Carina, Sculptor, and Leo II). The fact that these two 
lowest bins are not consiste nt with zero galaxie s has im- 
portant implications for the lStoehr et alj (|2002f ) solution 
to the MSP: specifically, it implies that the 11 well-known 
MW satellites do not reside in subhalos that resemble the 
11 most massive subhalos in Via Lactea. 

To further emphasize this point, we see from Figure [5] 
that the mass of the 11th most massive subhalo in Via 
Lactea is 4 x 1O 7 M . From the likelihood functions in 
Figure [TJ Sextans, Carina, Leo II, and Sculptor must 
have values of M . 6 less than 4 x 1O 7 M at 99% c.L, 
implying a negligible probability that all of these dSphs 
reside in halos with M .6 > 4 x 1O 7 M . 

Using the Mo. 6 mass function of MW satellites, we 
can test other CDM-based solutions to the MSP. Two 
mod els of interest are based on reioniza tion suppres- 
sion ([Bullock et al.ll2000tlMoore et al.ll2006ft and on there 
being a c haracteristic halo mas s scale prior to subhalo 
accretion ([Diemand et al.l feOOTf ) . To roughly represent 
these models, we focus on two subsamples of Via Lactea 
subhalos: the earliest forming (EF) halos, and the largest 
mass halos before they wer e accreted (LBA) int o the 
host halo. As described in iDiemand et al.l (|2007h , the 
LBA sample is defined to be the 10 subhalos that had 
the highest V^ lax value throughout their entire history. 
These systems all had V max > 37.3 km s -1 at some point 
in their history. The EF sample consists of the 10 sub- 
halos with Vmax > 1 6.2 kms" 1 (the l imit o f atomic cool- 
ing) at z = 9.6. The lKravtsov et all (|2004h model would 
correspond to a selection intermediate between EF and 
LBA. In Figure [5] we show the observed mass function of 
MW satellites (solid, squares) along with the EF (dotted, 
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Fig. 6. — The solid and dashed curves show the MW satellites 
and dark subhalos in Via Lactea, respectively. These lines are 
reproduced from Figure \5\ with error-bars suppressed for clarity. 
The blue (dotted) curve represents the ten earliest forming halos 
in Via Lactea, and the green (long-dashed) curve represents the 10 
most massive halos before accretion into the Milky Way halo. 

triangles) and LBA (long-dashed, circles) samples. We 
conclude that both of these models are in agreement with 
the MW satellite mass function. Future observations and 
quantification of the masses of the newly-discovered MW 
satellites will enable us to do precision tests of the viable 
MSP solutions. Additionally, once the capability to do 
numerical simulations of substructure in warm dark mat- 
ter models becomes a reality, the Mo. 6 mass function will 
provide an invaluable tool to place constraints on WDM 
models. 

5. SUMMARY AND DISCUSSION 

We have provided comprehensive dark matter mass 
constraints for the 9 well-studied dSph satellite galax- 
ies of the Milky Way and investigated CDM-based solu- 
tions for the missing satellite problem in light of these 
constraints. While subhalo V max values are the tradi- 
tional means by which theoretical predictions quantify 
substructure counts, this is not the most direct way 
to confront the observational constraints. Specifically, 
Knax is poorly constrained by stellar velocity dispersion 
measurements, and can only be estimated by adopting 
cosmology-dependent, theoretical extrapolations. We ar- 
gue the comparison between theory and observation is 
best made using the integrated mass within a fixed phys- 
ical radius comparable to the stellar extent of the known 
satellites, ~ .6 kp c. This approach is motivated by 
iStrigari et alj (|2007l ) who showed that the mass within 
two stellar King radii is best constrained by typical ve- 
locity dispersion data. 

Using Mo. 6 to represent the dark matter mass within a 
radius of 0.6 kpc, we computed Ma. 6 likelihood functions 
for the MW dSphs based on published velocity dispersion 
data. Our models allow for a wide range of underlying 
dark matter halo profile shapes and stellar velocity dis- 
persion profiles. With this broad allowance, we showed 
that the M .6 for most dwarf satellites is constrained to 
within - 30%. 



We derived the M .6 mass function of MW satellites 
(with error bars) and compared it to the same mass func- 
tion computed directly from the Via Lactea substruc- 
ture simulation. While the observed Mo. 6 mass function 
of luminous satellites is relatively flat, the comparable 
CDM subhalo mass function rises as ~ M^l. We rule 
out the hypothesis that all of the well-known Milky Way 
satellites strictly inhabit the most massive CDM subha- 
los. If luminosity does track current subhalo mass, this 
would only be possible if the subhalo population of the 
Milky Way were drastically different than that predicted 
in CDM. However, we show that other plausible CDM so- 
lutions are consistent with the observed mass function. 
Specifically, the earliest forming subhalos have a flat Afo.6 
mass function that is consistent with the satellite subhalo 
mass function. This would be expected if the population 
of bright dwarf spheroidals corresponds to the residual 
halo population that accreted a sig nificant mount of ga s 
before the epoch of reionization (|Bullock et al.l [2000) . 
We also tested the hypothesis that the present dwarf 
spheroidal population corresponds to the subhalos that 
were the most massiv e before they fell into the MW halo 
(jKravtsov et al.ll2004f ). This hypothesis is also consistent 
with the current data. 

In deriving the Mo. 6 mass function for this paper we 
have set aside the issue of the most- recently discovered 
population of MW dwarfs. We aim to return to this 
issue in later work, but it is worth speculating on the 
expected impact that these systems would have on our 
conclusions. If we had included the new systems, mak- 
ing ~ 20 satellites in all, would it be possible to place 
these systems in the ~ 20 most massive subhalos in Via 
Lactea? Given the probable mass ranges for the new 
dwarfs, we find that this is unlikely. We can get a rough 
estimate of their masses from their observed luminosi- 
ties. We start by considering the mass-to-light ratios 
of the known dSph population from figure [3] and from 
iMateol fl998). If we assume that the other dwarfs have 
similar M/L range, we can assign a mass range for each 
of them. In all cases, the new MW dwarfs are approxi- 
mately 1 to 2 orders of magnitude smaller in luminosity 
than the well-known dSph population. Using the central 
points for the known dSphs, we obtain Mq.q/L spanning 
the range from 3 — 230. Considering the width of the 
likelihoods, we can allow a slightly larger range, 2 — 350. 
If we place the new dwarfs in this latter range, the uncer- 
tainty in their masses is (2 — 350)LMq/ L Q . Even with 
this generous range we expect most of the new dwarfs 
have M .6 < 1O 7 M . 7 

The discovery of more members of the MW group, 
and the precise determination of the Mo. 6 mass func- 
tion, could bring the status of the remaining viable MSP 
solutions into sharper focus. These measurements would 
also provide important constraints on warm dark matter 
models or on the small scale power spectrum in CDM. 
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